function dS = STM_CR3BP(~,Y,u)

currentStates = Y(37:42);
%[] State at current integration step

x = currentStates(1);
y = currentStates(2);
z = currentStates(3);
vx = currentStates(4);
vy = currentStates(5);
vz = currentStates(6);
%[] Current States.

Umatrix = State2U_matrix(currentStates,u);
%[] Determine the U matrix. (Partial Partial derivatives)

K=[0, 2, 0;
  -2, 0, 0;
   0, 0, 0];
%[] Form K matrix

Df = [zeros(3),eye(3);...
      Umatrix   ,   K];
%[] Determine the 

A = zeros(6);
for i=1:6
    for j=1:6
        A(i,j)=Y(6*(i-1)+j);
    end
end
% Current STM matrix values.

DfA_temp=Df*A;
%[] Determine the time derivatives of the state transition matrix.

DFA = zeros(36,1);
%[] Pre-allocate memory for time derivative of the STM

for i=1:6
    for j=1:6
        DFA(6*(i-1)+j)=DfA_temp(i,j);
    end
end  
%[] change the _temp, which is the matrix to DFA, which is expressed in
%column vectors

Xdot = zeros(6,1);
%[] Preallocate memory

r1 = ((x+u)^2+y^2+z^2)^(1/2);
r2 = ((x+u-1)^2+y^2+z^2)^(1/2);
%[] Range of the satellite from first and second primary bodies.

Xdot(1:3) = [vx;vy;vz];
%[] Time derivative of the position is the velocity

Xdot(4)= x  -  (1-u)*(x+u)/r1^3  -  u*(x+u-1)/r2^3 + 2*vy;
Xdot(5) = y  -  (1-u)*y/r1^3 - u*y/r2^3 - 2*vx;
Xdot(6) = -(1-u)*z/r1^3-u*z/r2^3;
%[] acceleration expressed in the CRTBP frame.

dS = [DFA;Xdot];
%[] Output differentials\




end